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ABSTRACT 


Electromagnetic scattering from trees and vegetation is of prime importance in 
radar and remote sensing. The actual problem of scattering from trees is rather 
complicated and involves three dimensional scattering from lossy, electrically large, 
and randomly oriented objects. 

In this thesis, the radar cross section of a planar fractal tree is considered. 
Although a planar tree is far from being real, scattering from it sheds light on the 
scattering phenomenon from an actual tree. The planar tree is generated using 
fractal geometry and its branches are considered perfectly conducting. The tree is 
illuminated by a plane wave and the problem is solved using the moment method. 
Data is presented for the radar cross section for different branching angles of the 


tree and at different frequencies. 
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I. INTRODUCTION 


A. NEED FOR THE STUDY 

The trees existing in the natural world are fractal anisotropic. They are made 
up of long, intersecting and lossy objects. The geometry of these objects is not easy 
to set up as they are randomly oriented in a three dimensional space. The use of 
fractals facilitates the modeling of semi-—randomly distributed structures. 
Mandelbrodt [Ref. 1] has shown that a number of naturally occurring phenomenon 
such as coastlines, clouds, trees, etc. are fractal in nature. For instance, when a 
branch is divided into two (or more), the ratio of the length of the subbranches to 
the main branch length remains constant. Furthermore, the branching angle also 
remains the same. 

In this thesis, scattering from planar fractal trees is considered. The radar 
cross section of a real fractal tree is a complicated scattering problem. The analysis 
of this problem in a two dimensional space gives an approximate idea of the radar 
cross section of a real tree. 

The radar cross section of an object is a quantitative measure of the ratio of 
the power density that is received and scattered by the object to the power density 
of the electromagnetic wave that illuminates that object. The radar cross section is 
independent of the range of the object for the far—field situation. The theoretical 


definition of the radar cross section is given by the formula: 


P= aa lim | = 
R- ow 


where E! is the incident electric field vector, E° is the scattered electric field vector, 
and R is the distance between the scattered object and the point of observation. The 
radar cross section has dimensions of area. Usually, it is expressed in square 


wavelengths. 


B. STATEMENT OF THE PROBLEM 

This thesis investigates the radar cross section of planar fractal trees. These 
trees are composed of planar thin strip dipoles of arbitrarily dimensions and 
Orientations in a plane. These structures are excited by plane waves of various 
frequencies. 

The first step in solving the problem is to calculate the induced current 
distribution on each strip. The calculation of the current distribution is based on the 
theory of the moment method and requires a knowledge of the impedance between 
any two of these strips as well as the voltage on each planar strip due to the 
incident electric field. The basic concepts and the calculation of the current 
distribution are described in Chapter 2. 

In Chapter 3, the development of a FORTRAN program, is presented. The 
evaluation of the radar cross section requires the knowledge of the scattered electric 
field due to the induced currents on the planar strips. The program computes the 
scattered electric field and then the radar cross section of that structure. The details 
of these calculations are also presented in this Chapter. 

The computer models that are used by the developed program are presented in 
Chapter 4. Their generation is based on the fractal geometry. An existing and 
modified program is used to generate the geometry of the planar fractal trees in 


order to be used as input in the developed program. 


The numerical results of the radar cross section of a single planar dipole and a 
a number of planar fractal trees are presented in Chapter 5. The scattering from a 
single dipole is compared with standard results for a similar case. The limitations of 
the developed program are also presented in this Chapter. 

In Chapter 5, the conclusions of the radar cross section results and 
recommendations are presented. The two programs that are used for this 


investigation are listed in the Appendices. 


Il. METHOD OF MOMENTS THEORY 


In this section of the study, the basic concepts of the moment method theory 
are presented. This theory is used in the development of the RCS program to find 
the current distribution on a planar strip due to an incident plane wave. 

For a given structure consisting of planar dipoles the impedance between any 
two of them is calculated from the knowledge of the geometry and the wavelength of 
the incident plane wave. The voltage on each dipole is calculated from the 
knowledge of the characteristics of the incident electric field. The induced current 
distribution on each dipole of the structure is determined from the calculated 


impedance and voltage using the method of moments theory. 


A. GENERAL THEORY 
The method of moments is a numerical procedure for solving integral 


equations of the form: 
b 
| iGc ) Geax \dx 7 ex) edie 5 (eqn 2.1) 
a 


where f(x’) is an unknown function, K(x,x’) is a known Kernel or Green's function, 
and g(x) is a given function. This procedure reduces the integral equation (eqn 2.1) 
to a system of simultaneous linear algebraic equations in terms of some unknown 
coefficients. This method requires that the function f(x’) be approximated by a 


series of N expansion functions or " basis functions ", such that 
3 


eS es SE anfalx’) Wo", yee nN (eqn 2.2) 
where the domain of f,(x’) is the same as that of f(x’) and ay's are the complex 
unknown expansion coefficients. 

There are two types of basis functions. The subdomain functions, which are 
nonzero over a part of the domain of the unknown function f(x’), and entire domain 
functions being nonzero over the entire domain of f(x’). In antennas, some 
commonly emploved subdomain basis functions are the piecewise sinusoid functions, 
the unit height—pulse functions and the piecewise triangular functions. 

The subdomain procedure requires subdivision of the structure into N 
nonoverlapping segments. Figure 2.1 shows a segmented line where the segments are 


assumed to be collinear and of equal length, although this condition is not necessary. 


SUBINTERV AL 


Figure 2.1 Segmented Line. [From Ref. 2] 


CA 


Figure 2.2 shows a subdomain unit height—pulse function which produces a 
staircase representation of the unknown function f(x’). Figure 2.3 shows a 
subdomain sinusoid basis function and the representation of the function f(x’). 
Figure 2.4 shows a subdomain triangular basis function producing a smoother 
representation of the function f(x’) than the case of the unit height—pulse basis 
function. 

The use of entire—domain basis functions does not require any segmentation of 
the structure. One of the most most commonly used basis functions of this kind is 


the sinusoidal basis functions. 
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Figure 2.2 Unit Height—Pulse Basis Function and a Staircase 


Representation of f(x’). [From Ref. 2] 





Figure 2.3 Sinusoid Basis Function and a Represeutation 


of f(x’). [From Ref. 3] 
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Figure 2.4 Triangular Basis Function and a Representation 


of f(x’). [From Ref. 2] 


The substitution of the function f(x’) by a sequence of N basis functions leads 
to one equation of N unknowns of the expansion coefficients an, which can be found 
by using N linearly independent equations. These equations are set up by the use of 
"testing or weighting" functions w»(x), m=1,2,....N. At this point the definition of 
the inner product is required. The inner product (f, g) between any two functions or 
vectors f, g is a scalar operation defined as 

Gg) = ie f-g ds (eqn 2.3) 
where S is the surface of the structure that is analyzed [Ref 4]. The inner product of 
the selected testing functions w(x), with the two sides of the original integral 


equation leads to the equation: 


( wal), i f(x’ )K(x,x" dx’) = ( g(x), wn() ) 


(eqn 2.4) 


or 


(ato 


The use of the above definition for the inner product, yields: 


N 
Y anfa(x’ )dx’ )= ( ax), W(x) ) eee 
=1 


n 


N b b b 
y an| un(x)ax | fae (ane doe! = | Scan Ox. ile lee 
=1 a a a 


n 


10 


The following substitutions 
b 
y= | g(x) wn(x)dx, (00 alee N (eqn 2.7) 
a 


and 


b b 
Ze | un(x)ax | foe eee) deems =1,2,....N 
a a 


(eqn 2.8) 
result in a matrix equation of the form: 
[Zmn] [aa] = [Vo] (eqn 2.9) 


The unknowns [ap] can be obtained through a matrix inversion: 
-1 
[an] = [Zma] - [Vn] (eqn 2.10) 


where Zn) is the inverse matrix. The column vector [V,] depends upon the given 
function g(x) and the selected testing functions w(x). The matrix [Zn] depends 
upon the known kernel K(x,x’) and both the selected basis and testing functions. 
Once the expansion coefficients are known, the function f(x’) is also known. 

The choice of basis and testing or weighting functions is based upon experience 
and the rule is that their number has to be the same. The procedure of using the 


same basis and testing functions is called the " Garlekin's method ". 


B. APPLICATIONS TO EM THEORY 

A number of problems in electromagnetic radiation and scattering can be 
solved by the method of moments. In the present case, the simple case of a perfectly 
conducting object " situated in free space " is considered. 

The basic problem is to investigate the case when an object is illuminated by 


fields of known impressed electric and magnetic currents Ge M)), In the absence of 


bE 


the object, the impressed currents radiate the assumed known incident electric and 
magnetic fields (EL, 1). In the presence of this object, the impressed currents 
radiate the unknown total fields (Ee H’). 
The integral equation is obtained by the surface equivalence principle, 
replacing the object by free space together with the electric surface current density 
J=nxit (eqn 2.11) 
where J exists on the entire surface S of the object and n is the unit vector normal 
to that surface. In free space, J radiates the scattered fields ES, H 
Ee (eqn 2.12) 
nH = Ht (eqn 2.13) 

The boundary conditions for this case enforce the total tangential electric field 

on the surface S to zero: 

nx(E°+E)=0 (eqn 2.14) 
This is an integral equation for J since the sc:ittered electric field E® can be written 
as an integral over S of the dot product of J and the dyadic free space Green's 
function. [Ref. 5] 

As the geometry of the object is known, it is more convenient to use the total 
current I instead of the total current density J. So, the problem is to find the 
unknown current I induced by the incident electric field. The method of moments 
solve these kind of problems by the following procedure. 

The first step is to expand the unknown current I in terms of some basis set: 

N 
ts eben (eqn 2.15) 


n=1 
where the I, are the sequence of N unknown complex coefficients, and the Fy is a 


12 


sequence of N known modes or basis functions. The best choice of F, for a given 
problem could be quite involved and is discussed in (Ref. 4, pp 308-310]. 

The second step is to select the testing or weighting functions w,, m=1,2.....N. 
These can be identical with the basis functions or different. The inner product of the 
sequence of N weighting functions wz, with both sides of the integral equation gives 
a N x N system of simultaneous linear algebraic equations of the symbolic form: 

(Z]-{1] = [V] (eqn 2.16) 

where I is the current column vector whose N components give the values of I, 


and [Z] is the N x N impedance matrix given by the equation 
Tee =— -|| ES. wpds (eqn Dal) 
5 


where S is the surface of the structure being analyzed. The impedance matrix [Z] is 
always symmetric, and, for the special case of a thin dipole instead of an arbitrary 
object, is also a Toeplitz matrix. In a Toeplitz matrix, Zm_ depends only on |m-—n|. 

Generally, [Z] is dependent only on the geometry and material composition of 
the scatterer, but not on the incident fields [Ref 5]. The right—hand side of the last 
equation is the voltage vector whose N components give the corresponding mode 
voltage. The voltage vector depends only on the excitation, i.e., the incident electric 
field. The dimensions of the elements of {Z] and [V] are volt-amps (VA), while the 
elements of {I] are dimensionless. 

The solution of the last matrix equation is the column vector [I], whose 
elements represent the complex coefficients Ip. As the total current I was expanded 
by a sequence of N known modes or basis functions and the complex coefficients Ip 


are known, the total current and so the total current density J is also known. 
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Although the choice of weighting functions is free, it has to be considered that 
the matrix equation being solved requires the evaluation of N? terms and each term 
requires two or more integrations. When these integrations are to be done 
numerically, the computations become complicated. There is a way to reduce this 
complexity by choosing as weighting functions the Dirac delta functions. This is the 
method of point—maiching in which delta functions are enforced only at discrete 
points on the surface S. The results of this method can be quite accurate especially 
when the discrete points are selected to be equally spaced. The solution satisfies the 
electromagnetic boundary conditions (e.g., vanishing tangential electric fields on the 
surface of an electric conductor) only at discrete points [Ref 4]. Between these 
points the boundary conditions may not be satisfied. In this case it is required to 
define the deviation as a residual (e.g., residual=AE|tan = Bere E' tan # 0 on 
the surface S of an electric conductor) and use the method of weighted residuals so 
that the boundary conditions will be satisfied in an average sense over the entire 


surface S. 
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I. ANALYSIS AND DEVELOPMENT OF RCS PROGRAM 


In this chapter, the basic steps that are involved in the calculation of the 
Radar Cross Section (RCS) of a planar structure composed of conducting strips are 
presented. As mentioned earlier, the fractal tree is modeled as consisting of planar 
thin strips. The radar cross section of the fractal tree is calculated using the method 
of moments. A Fortran program is developed to calculate the RCS of the tree for a 
specified geometry and at a given frequency. 

The moment method discretizes an integral equation to a matrix equation of 
the form: 

(Z}- {1} = [Vj (eqn 3.1) 
where [Z] is the impedance matrix whose elements represent the mutual impedance 
between any two dipoles of a given structure depending upon the wavelength and 
geometry of the structure, [V] is the voltage matrix whose elements correspond to 
the voltage on each dipole due to excitation ( incident electric field ), and [I] is the 


unknown matrix whose elements represent the induced current on each dipole. 


A. DEVELOPMENT OF RCS PROGRAM 

The numerical results for RCS are obtained from a Fortran program that 
calculates backscattered RCS of a planar structure consisting of arbitrarily oriented 
dipoles. For convenience, the structure will be assumed to be in the y—z plane of an 
xyz cartesian coordinate system. 

Figure 3.1 shows the structure to be investigated. A large number of planar 


dipoles of variable lengths, widths, and orientations in the y—z plane is 
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illuminated by a plane wave with electric field linearly polarized and characterized 
by the angles yo and @p. It is required to calculate the backscattered RCS from this 


structure. 


The incident electric field E' is given by the formula 


Bi = e.g Hkx-x + ky-y + ke-z) (eqn 3.2) 
where k is the propagation vector and 
kyo + eae k." = ko" = Ww toed (eqn 3.3) 
ky = el (eqn 3.4) 
r 
kx = kosin@cos Yo (eqn 3.5) 
ky = kosin@psin yo (eqn 3.6) 
k, = kocos@ (eqn 3.7) 
= Exx + Eyy + Ez (eqn 3.8) 


The electric field Bi lies on the plane perpendicular to the direction of 
propagation of the plane wave. Hence, eal = 0 implies that 

Exkx + Eyky + Ezk, = 0 (eqn 3.9) 
The substitution of Ey and E, by the variables a and b (independent variables) gives 


the coordinate Ex as a dependent variable 


—— (aky + bk) (eqn 3.10) 


kx 
Equation 3.1 the becomes 


kx 
(eqn 3.11) 
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Figure 3.1 Geometry of RCS Program. 


def 


The RCS program that has been developed makes the following assumptions: 
1. The dipoles are very thin planar strips. The width of the dipoles is 
assumed to be electrically small so that only the axial current is 
significant. Further, there is no current variation along the width of 
the dipoles. 
2. The dipoles are not intersecting. 
3. The dipoles are perfect conductors. 

The program uses the theory of moment method and the way that it solves the 
problem can be understood if a single planar thin dipole in the y—z plane is 
considered. Each dipole is subdivided into equal segments. Overlapping Piecewise 
Sinusoidal (PWS) modes are assumed to exist on the dipole. The length of each 
PWS mode is equal to the length of two segments. Figure 3.2 shows the mth PWS 
mode. The PWS mode has a length 2h,, and a width 2w». The coordinates of its 
center are (ym, Zm), and it is oriented at an angle yi measured from the z axis. 

The incident electric field induces current along the axis ¢ of that mode. In 
Figure 3.2, the induced current J varies sinusoidally along the axis of the my, mode. 
A new coordinate system 7~-¢ is introduced, where ¢ is the axis of the dipole and 7 is 
the axis perpendicular to ¢. The 7-¢ coordinate system is obtained by rotating the 
y—z system about the x—axis through an angle 90°—y,. Figure 3.3 shows the details 
of this rotation. 

The relation between the coordinates of the center of the mode along the axis 
of the two orthogonal systems is found to be: 

y = ym + Csin( Ym) — nC08( Ym) (eqn 3.12) 
Z = Zm + Ccos(Ym) + nsin( Ym) (eqn 3.13) 


18 


The corresponding vector equation is: 
y = ¢sin( Um) — 7cos( vn) (eqn 3.14) 


z = ¢cos(vm) + 7sin( Yn) (eqn 3.15) 





Figure 3.2 Geometry of the mzh PWS Mode of the Structure. 
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Figure 3.3 Transformation of PWS Mode's Center Coordinates. 


Ie Voltage Equations 
In this section, expressions for the elements of the voltage matrix [V] 
due to the incident electric field, are presented. Each element of this matrix 
corresponds to a PWS mode of the structure that is investigated. The size of the 
voltage matrix is equal to the total number of modes of the structure. 
The induced current density on the mth mode of the structure, due to 


incident electric field, is 


where 2wy is the width and 2h, is the length of the myn mode. For each mode, the 
induced current variation will be sinusoidal along its axis, being maximum at the 


center and zero at the end points. The corresponding voltage Vy is 


Wm Dm; 
V = | Bi. Jaan (eqn 3.17) 
mm —Wp —hp xX 


From equations 3.11 and 3.16, it is seen that 


= [a sin(Ym) + b cos(vn)] eilky(ym + ¢sin(Yn) + ke(zm +6cos(Ym)] 


(eqn 3.18) 
where Wp is the angle between the z axis (reference for angle measurements) and the 


axis of the mode and yn, 2m are the coordinates of the center of the mode along the 


al 


y and z axis respectively. The substitution into equation 3.17, gives the following 


form of Vn: 
ee ekyym + kezm) pad 
Ma = fast) +b on) |, HS sin(kolte | ¢D)a¢ 


(eqn 3.19) 
where ke = kysin(~m) + kzcos(¥%m). A closed form evaluation of the integral in this 


equation leads to final form of Vn 


Vn = —pkoVon— | cos(kobin) — cos(k cha) ke ## ko 
ke re ko 

(eqn 3.20) 

Van = Von hn sin(kohm) ke = # ko 
(eqn 3.21) 

where 
—j(kyym + kz2m) 
V.  ———r | a sin(um) + b cos( qm) 
sin(kohn) 

(eqn 3.22) 

2. Radiation Equations 


In this section, expressions for the far—zone scattered fields due to the 
mth PWS mode are developed. For far—field observations the electric field is given 


in spherical coordinates by the following equations [Ref 4]: 
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‘* - (eqn 3.23) 
s]- —jkor 


4 ? @) (eqn 3.24) 
aT 
ike -jkor 
Beet (kg aN.) (eqn 3.25) 
rey 


where 


Ly= N [Mxcos@cosy + Mycosésiny —Mzsin 6} et dkor’costm go/ 


(eqn 3.26) 
L. a I [-Mxsiny a Mycosy et dkor “cos Yn ds’ 
S 


(eqn 3.27) 
Ng = {| [Jxcos@cosy + Jycosésiny — Jzsin}] et Jkor’ cos} n fa 
5 


(eqn 3.28) 
Me = | [—Jxsiny 7 Jycos 7} et Jkor COS Wm Ae 
S 
(eqn 3.29) 
7 — 1205 (eqn 3.30) 


The quantities Jx, Jy, Jz, are the components of the electric current 
density De that are induced on the mt, mode over the surface S, and the quantities 
Mx, My, Mz, are the coordinates of the magnetic current density M, over the surface 
S. As the structure is on the y—z plane and Mg is zero, the quantities Jy, L 9 L Y are 


zero. Equations 3.24 and 3.25 become: 


‘Le JKor 
E x tikoe “9 N, er 
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As the current density deyis along the ¢—axis, all quantities within the 
integrals for J and M will be expressed in terms of the 7—-¢ system. The surface 
element that is used in all integrals is ds’ = dydz = d¢dn = d¢, as the width of the 
mode is assumed to be very small in terms of its length. The transformations from 


the y—z system to the (—7 system are made by using the following substitutions: 


t’cos(Um) = ysin&iny + zcosé (eqn 3.33) 
where 

Y= Yn + Csin(yn) (eqn 3.34) 

Z = Zm + Cc08( Yn) (eqn 3.35) 
and 

In = Jy + Juz = ¢ | —42 — sin(bo(he-1¢1))) 

2wnsin(kohm) 
(eqn 3.36) 

where 

Jy = Je sin( tm) (eqn 3.37, 

Jz= Je cos( tm) (eqn 3.38) 


These substitutions and algebraic manipulations give the quantities Np, and N om 


for the m¢h mode in the (—n system: 


N gin = Ng 32g | €08(Enbn) — 008(Kohn) En? 2ko 
Bn = ko 
(eqn 3.39) 
Nom = Ng JG sin(kghp) En = +ko 


(eqn 3.40) 


and 
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, 1 2k 
ea ="N ms | cos(Emhm) — cos(kohn) En? to 
Ekg 
(eqn 3.41) 
Nom — Ms jie sin(Enh;) ce = +ko 
(eqn 3.42) 
where 
Ny = sinvimcos@siny — cosemsind a ollko(ynsin &iny+zmcos 8)] 
sin(kohn) 
(eqn 3.43) 
N = Sinvncosy | ellko(ynsin &iny+2mcos8)| 
ee mM 
YO sin(kolin) 
(eqn 3.44) 


En = Ko(sinypsin &iny + cosypcos 4) 
(eqn 3.45) 


In this thesis, only the monostatic radar cross section will be 


considered. In the radiation equations the angles @ and y represent the orientation 


angles of the scattered electric field E°. Their relation with the incident angles 


and ¥ for the monostatic case is: 
6=7—% (eqn 3.46) 
Y= TY (eqn 3.47) 


If M denotes the total number of PWS modes in the structure under 


investigation, then 
(eqn 3.48) 


M 
Nae oe ON 


(eqn 3.49) 
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The final expression for the equations 3.24 and 3.25 is 


M 
E,—c: Sen (eqn 3.50) 
8 ia 6m 
¥ 
a N , 
y —— ym (eqn 3.51) 
where 
ap oJ kor 
C= = aaa (eqn 3.52) 
Arr 
3. RCS Equations 


When an incident plane wave with electric field 5’ strikes the object 


and F° is the scattered electric field, the radar cross section o is defined as 


=S 
o = lim 41” El. (eqn 3.53) 


For the scattered electric field e its spherical coordinates E ° and a have been 
calculated by the equations 3.50 and 3.51. The incident electric field E’ is known by 


means of equation 3.11. 





ee 2 2 2 
[sie alaiiee al ieee St (eqn 3.54) 
Kx 
oie 2 2 2 2 2 
= E = 
I] =1Egl + 1B =ICl [LINgl +INyI | 
(eqn 3.55) 
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The combination of equations 3.11, 3.50—3.52, 3.54—3.55, and 3.53 leads to the final 
formula for RCS: 
2 I 


2 2 Dy 
ba RB) + aky + bk, _ 


x 








Baie 
o = S00 | |Ngl? + [N 
An G y 


(eqn 3.56) 


This is the formula that the program uses to compute the RCS of a 


planar structure for a given set of incident angles 6, and Yo. 


B. INPUT-OUTPUT OF RCS PROGRAM 
The RCS program, which is described in Appendix A, is a FORTRAN 
program having two input files. The first input file, INPUT1/, contains the data that 
characterize the incident plane wave and the data that describe the geometry of the 
planar structure whose broadside RCS is measured. The second input file, INPUT, 
contains the set of incident angles 6, Yo. 
The program reads from the INPUT} file the following input data: 
1. frequency of the incident plane wave in GHz, 
2. parameters a and b that characterize the polarization of the incident 
electric field. 
number of dipoles that the structure consists of, 
half length of each dipole in cm, 
half width of each dipole in cm, 


Cee oe 


coordinates in cm of the center of each dipole along the two axes of 


the orthogonal system that the structure lies, 


ra 


7. orientation angle of each dipole, measured from the vertical axis, 
positive in the clockwise direction, and 
8. number of segments of each dipole. 

As the program reads these input data it generates the geometry of the PWS 
modes, calculates the impedance elements between any two modes of the structure, 
and fills the impedance matrix [Z]. The size of the impedance matrix depends on the 
total number of modes. When this matrix is calculated and filled, the program 
inverts it to Zi and stores it for later use. 

The next step is to read from the INPUT file the sets of incident angles 4%, Yo 
and for each one it calculates the voltage on each mode, filling the voltage matrix 
[V]. This matrix is a column vector which depends upon the excitation only. Then it 
multiplies the stored inverted impedance matrix in by the voltage matrix. This 
yields the current vector [I]: 

] = (z—1-{v] (eqn 3.57) 

As the induced currents are known, the program uses the previously described 
radiation equations and calculates the RCS of the structure corresponding to the 
given set of incident angles 6, Y. The output RCS is normalized to the square of 
the wavelength ae and is given in dB. The program reads the next set of incident 
angles 6, Y and repeats the same procedure to compute the RCS of the new set. 

Although the input lengths and widths are in cm, the program considers them 
normalized to the wavelength. For accurate results at least 4 segments per 
wavelength are chosen. The selection of the width of each dipole is arbitrarily taken 


as L/W = 33. 
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IV. FRACTAL TREE GENERATION 


The problem of measuring the radar cross section of a real natural tree is 
complicated as it requires an investigation in three dimensional space with very 
long, intersecting elongated objects that are randomly oriented. 

For simplicity the radar cross section (RCS) problem is solved in a two 
dimensional space. RCS is calculated from a planar fractal tree whose branches are 
considered as thin planar dipoles of different lengths and widths. In an actual tree, 
the branches are lossy and in general anisotropic. However, in the present model, 
the tree is comprised of perfectly conducting branches. The geometry of this tree is 


based upon the fractal geometry. 


A. DEFINITIONS 

Fractal is a mathematical set or object whose form is extremely irregular or 
fragmented at all scales [Ref. 6]. The requirement to describe the shape of many 
objects that appear in the natural world, such as trees, mountains, coastlines, etc., 
led to the generation of fractal geometry. As Euclidean geometry cannot give 
mathematical expressions to describe fractal objects, fractal geometry is used to 


describe mathematically many natural patterns. 


29 


The basic characteristics of fractal objects are [Ref. 7]: 

1. A large degree of heterogeneity. 

2. <A self-similar structure over many size scales. Self—similarity refers 
to the general preservation of form or characteristic regardless of the 
scale of observation. 

3. The lack of a well-defined (characteristic) scale. 

The geometric characteristics of fractal objects are useful for describing phenomena 
of nature, such as scattering of objects like landscapes or surface cracks. Although 
these objects are irregular and randomly oriented in nature, they show structural 
similarities on several different discrete size scales [Ref. 7]. 

One measure of structural complexity is the fractal dimension Dr. There are 
several definitions of Dr depending upon the particular application. For the case of 


fractal trees the fractal dimension Dr is the measure of the space that a 


self—similar structure filis, and it varies with the branching levels that the structure 
consists of. 
The most useful terms from fractal geometry that are used to describe a planar 
fractal tree are the following: 
1. The number of branch segments N formed from each preceding branch 
segment. 
2. The constant similarity ratio "r" that relates the fractional reduction 
in segment length for each segment to previous level. this factor is less 
than 1.0. 


In this case, the fractal dimension Dr is given by the formula [Ref. 7]: 


1—Dr = ogttotal branch length) _ _log(rN) 


log(average branch length) log(1/r) 
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This equation is accurate for symmetric structures or asymmetric structures having 
a normal distribution function [Ref. 7]. More details about fractal geometry are 


discussed in [Ref. 1] and [Ref. 7]. 


B. FRACTAL TREE GENERATION 

In the planar fractal structures that are used for calculating the radar cross 
section, the number of branch segments N formed from each preceding branch 
segment is 2. The principal properties of these trees are that the branches are not 
overlapping, and the angle # between any two branches is the same. The 
nonoverlapping between the branches of each fractal model is achieved by choosing 
the branch angle @ for a given reduction factor by the following empirical formula 
[Ref. 7] 


irene = 3294 po! 


(Radians) 

This angle is given in radians and r is the desired reduction factor which is a 
constant for the structure. As the reduction factor r increases the tree fills more 
space. 

Five types of planar fractal trees are considered in this thesis. Each type 
corresponds to a set of values for reduction factor r and minimum branch angle @. 
Table 4.1 shows these sets of r and minimum 8, where @ is given in degrees. 

Each model is generated by selecting the desirable reduction factor r (r < 1.0), 
the corresponding branch angle @ given by equation 4.1, the initial length from 


which the reduction will start, the value of N (which in the investigated case is 2), 
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and the number of the desired levels for branch generation. Each level contains 2° 


points corresponding to the end points of the reduced length branches. 


TABLE 4.1 
NUMERICAL VALUES OF r AND MINIMUM ANGLE @. 


REDUCTION FACTOR (1) ANGLE @ 
0.53 29.940 
0.55 46.430 
0.60 95.89° 
0.66 145.350 
0.71 180.00° 


Figures 4.14.5 show the planar fractal trees that correspond to each set. of r 
and @ of Table 4.1 respectively. As the branch angle @ increases, the spreading of the 
branches becomes larger. In all trees the physical length of the initial branch is 
chosen as two centimeters. These models were generated by a FORTRAN program 
which has the following input data: 

1. The initial length. This the length of the first dipole whose length will 
be reduced by the constant reduction factor (r). 

2. The value of N, being 2 in all cases. 

3. The initial point that the reduction starts. This point is the end point 


of the initial wire. 
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4. The number of desirable levels for branch generation. 
5. | The branch angle @. 

The program generates two output files. The first is the input file for the RCS 
program containing all the required data that were mentioned in Chapter 3. The 
second output file contains the coordinates of the start and end points of each 
generated branch. 

For the models characterized by reduction factors 0.53, 0.55, and 0.60, the 
program limits the structure size when the length of a branch becomes smaller than 
0.1 of the wavelength of the incident plane wave. For the models characterized by 
reduction factor 0.66 and 0.71, this limit is taken as 0.15. The reason is that as the 
reduction factor increases the number of branches in a particular branch level 
increases. The size of the tree is truncated so that the number of unknowns is 
manageable. 

The geometry of all planar fractal trees is in the same coordinate system that 
the RCS program uses. All models have been generated in the y—z plane of a 


cartesian coordinate system. 
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Figure 4.1 Fractal Tree for r = 0.53 and @ = 29.94°. 
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Figure 4.2 Fractal Tree for r = 0.55 and @ = 46.439. 
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Figure 4.3 Fractal Tree for r = 0.60 and @ = 95.89°. 


36 


Figure 4.4 Fractal Tree for r = 0.66 and @ = 145.35°. 
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Figure 4.5 Fractal Tree for r = 0.71 and @ = 180°. 
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V. NUMERICAL RESULTS 


In this chapter, computed results for the radar cross section o of planar 
structures are presented. In order to check the program, the radar cross section 
results of a single planar dipole are compared with the results given in [Ref. 8]. 
Details of this comparison are presented in the following section. 

In all models where numerical results for the radar cross section are presented, 
the following factors have been considered: 

1. E-—plane is the plane corresponding to yp = 0°, and 6 varying from 0° to 
180°. 

2. H—plane is the plane corresponding to % = 90°, and yo varying from —90° 
to +909. 

3. The resulting RCS o is normalized to the square of the wavelength A of the 
incident plane wave (0/2). 

4. The number of segments, that each dipole is subdivided, is taken as four 
per wavelength. 

5. The physical dimensions and orientations of the dipoles used to construct a 
planar structure are the same when this structure is investigated at different 
frequencies, and only the segmentation is different depending upon the wavelength A 
of the incident plane wave. 

6. Each PWS mode has length equal to the length of two segments. 

7. The incident electric field is linearly polarized and the parameters a and b 
are taken as 0 and 1 respectively for this investigation. This corresponds to having 


only a horizontal magnetic field. 
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A. SCATTERING FROM A SINGLE DIPOLE 

In Figure 5.1 a centered loaded vertical planar thin dipole of length L and 
width 2w is oriented along the z axis. This dipole is excited by a plane wave of 
frequency 30 GHz traveling in the x—y plane (& = 90°, yo = 0°). The same 
situation is described in [Ref. 8, pp. 510-515] for a cylindrical dipole of radius a and 
length L such that L/2a=74.2. In the present case of a planar dipole, the length of 
the planar dipole is selected such that L/2w = 33 [Ref. 9]. Figure 5.2 shows, on a 
semi—log scale, the results computed by the developed program of the monostatic 
normalized radar cross section ep of this single dipole for values of L/A ranging 
from 0 to 1.4, and for loads ZL = 0 and ZL = w. The case of ZL = » is achieved by 
setting a small gap (0.01 wavelength) at the center of this planar dipole. In Figure 
5.2, the corresponding values of o/ d? from (Ref. 8, pp. 115] are also shown. 

This comparison leads to an accurate check of the correct calculation of the 
radar cross section by the developed RCS program. The small discrepanc:’ between 
the computed values and those given by [Ref. 8] is due to the error incurred in 


reading values off the curves presented in (Ref. 8, pp. 115}. 
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Figure 5.1 Vertical Centered Loaded Planar Dipole. 
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Figure 5.2 Normalized RCS o/ d” of a Center Loaded Planar Dipole. 
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B. SCATTERING FROM PLANAR FRACTAL TREES 

In Chapter 4, five types of planar fractal trees were described. Each type of 
these trees was lying in the y—z plane. The developed RCS program calculates the 
broadside radar cross section o of each tree for frequencies*15 GHz—75 GHz, for the 
monostatic case. The initial dipole that is being reduced by the reduction factor r 
has physical length 2 cm and only the branch lengths are changed for each tree 
depending upon the reduction factor r. As the frequency of the plane wave changes, 
the electrical length of this dipole as well as the dipoles composing the tree will also 
change. For the frequency range of 15 GHz-75 GHz, the electrical length of the 
initial length will vary from one to four wavelengths. 

Figures 5.3-5.6 show the variations of Epa in dB, in the E—plane and the 
H—plane for the case of a fractal tree characterized by reduction factor r = 0.53 and 
branch angle @ = 29.949. The frequency of the incident plane wave varied from 15 
GHz to 60 GHz in steps of 15 GHz. The tree is composed of 31 dipoles. This number 
is small as the reduction factor is small and the branch lengths are reduced to very 
small values at low levels. 

The of d variations in the E—Plane are in a range of 10 dB approximately. In 
the H—Plane. The o/ d? is symmetric about the 90° axis and is smoother at lower 
frequencies. As the frequency increases, the maximum value of o/ d? at § = 90° and 
Yo=0°, increases from 2.73 dB at 15 GHz to 19.23 dB at 60 GHz. Figure 5.7 shows 
the variation of the maximum Be in terms of frequency increments for the same 


structure. This variation is very small between 30 and 45 GHz. 
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Figure 5.3 Normalized RCS in dB at F = 15 GHz of a Fractal Tree 
with r = 0.53 and @ = 29.940. 
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Figure 5.4 Normalized RCS at F = 30 GHz of a Fractal Tree 
with r = 0.53 and @ = 29.940. 
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Figure 5.5 Normalized RCS at F = 45 GIlz of a Fractal Tree 
with r = 0.53 and @ = 29.949. 
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Figure 5.6 Normalized RCS at F = 60 GHz of a Planar Fractal Tree 
with r = 0.53 and #@ = 29.949. 


47 


i 
~m O) 
lO O) 
ON 
yd 
= (D 
j 
Ke} 
NR 
N 
So + 
OU 
+ > 
©) 
We 
Qi 
Oe 
© 
LW 
Sie 
ae 01 
@ 
wo @ LO eC Ww) @ 
WN WN Tea Na 
(ap) N/ Oo 
G 


Figure 5.7 Variation of maximum o/ d? vs. Frequency 


of a Planar Fractal Tree with r = 0.53 and @ =29.94°. 
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Figures 5.8—5.12 show the variations of ae , in dB, in the E and H planes 
for the fractal tree characterized by a reduction factor r = 0.55 and branch angle @ 
= 46.43°. In this case, the fractal tree is composed of 63 dipoles and spreads more 
than the previous one, and the lengths of the branches are bigger than the previous 
ones (due to a larger reduction factor of 0.55 instead of 0.53). This results in more 
variations for o/ \” in both E and H planes. 

This model is investigated at a frequency range of 15-75 GHz. At all 
frequencies, the o/ varies in a range of 10-12 dB approximately. In the H—Plane 
the variation of o/ ue is symmetric about the 90° axis. The maximum value of radar 
cross section varies from 2.78 db at 15 GHz to 22.63 dB at 75 GHz. Figure 5.13 
shows the maximum value of o/ 2 corresponding to # = 90° and yo = 0° in terms 
of the frequency of the incident plane wave, varying from 15 GHz to 75 GHz. This 
maximum o/ i increases linearly as the frequency increases except the range of 


45-60 GHz where the variation is very small. 
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Figure 5.8 Normalized RCS at F = 15 GHz of a Planar Fractal 
with r = 0.55 and @ = 46.43°. 
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Figure 5.9 Normalized RCS at F = 30 GHz of a Planar Fractal Tree 
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Figure 5.10 Normalized RCS at F = 45 GHz of a Planar Fractal Tree 
with r = 0.55 and @ = 46.439. 
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Figure 5.11 Normalized RCS at F = 60 GHz of a Planar Fractal Tree 
with r = 0.55 and 0 = 46.439. 
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Figure 5.12 Normalized RCS at F = 75 GHz of a Planar Fractal Tree 
with r = 0.55 and @ = 46.439. 
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Figure 5.13 Variation of Maximum RCS vs. Frequency 


of a Planar Fractal Tree with r = 0.55 and 6 = 46.439. 
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Figures 5.14—5.18 show the variations of the normalized radar cross section 
a/ re of a planar fractal tree characterized by reduction factor r = 0.60 and branch 
angle @ = 95.89°. This model is composed of 63 dipoles and has more spreading than 
the two models that were previously investigated. The branches have larger physical 
lengths than the ones in other two models. At 15 GHz, the variation of GAC in the 
E—Plane is very low and similar to the variation in the other two models at the 
same frequency. In the H—Plane this variation is different as the incident plane 
wave strikes more dipoles from —90° to 90°. 

As the frequency increases, the electrical length of the dipoles is also increased, 
and more lobes appear at 30 GHz and higher frequencies in both the E and the H 
planes. The maximum value of ape at = 90° and yo = 0° varies from —0.23 dB 
at 15 GHz to 21.60 dB at 75 GHz. Figure 5.19 shows the variation of the maximum 
radar cross section in terms of the frequency of the incident plane wave. In a range 
of 15-60 GHz the values of maximum o/ d? jollow a straight line approximately. In 
the range of 60—75 GHz the variation of maximum o/ d is small. 

The planar fractal trees characterized by r = 0.66, 6 = 145.359, and r = 0.71, 0 
= 180°, have large values of branching angles. The large values of reduction factor r 
generate fractal trees whose branches have large physical lengths compared with the 
lengths of the branches of the previously investigated fractal trees. The result is that 
the electrical lengths of these branches are large also at the range of 15—75 GHz, and 
a large number of modes is required to investigate the last two types of fractal trees. 
It was found that the developed RCS program is not able to calculate the radar 
cross section of fractal trees that are characterized by large values of reduction 
factor r and branching angle 6 due to memory restrictions and other numerical 


problems. 
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For this reason the numerical results for ax are not presented in this thesis for 
these two cases. 

It is seen by comparing Figures 5.7, 5.13, and 5.19 that there is a frequency 
range over which the variation of maximum RCS is rather small. Furthermore, this 
range of frequencies shifts to higher frequencies as the tree structure is spread from r 
= 0.53 to r = 0.60. The maximum RCS of each of these trees varies in a range of 3 
dB approximately at the same frequency. 

Scattering from an actual tree will not exactly follow all these patterns, but it 
is felt that the trends would generally remain the same. This is specially true for the 


frequency behavior. 
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Figure 5.14 Normalized RCS at F = 15 GHz of a Planar Fractal Tree 
with r = 0.60 and # = 95.899. 
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Figure 5.15 Normalized RCS at F = 30 GHz of a Planar Fractal Tree 


with r = 0.60 and @ = 95.899. 
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Figure 5.16 Normalized RCS at F = 45 GHz of a Planar Fractal Tree 
with r = 0.60 and @ = 95.899. 
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Figure 5.17 Normalized RCS at F = 60 GHz of a Planar Fractal Tree 
with r = 0.60 and @ = 95.899. 
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Figure 5.18 Normalized RCS at F = 75 GHZ of a Planar Fractal Tree 
with r = 0.60 and @ = 95.899. 
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Figure 5.19 Variation of Maximum RCS vs. Frequency 


of a Planar Fractal Tree with r = 0.60 and @ = 95.899. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The radar cross section of natural trees is a complicated problem of 
electromagnetic scattering. The real trees are made up of long, intersecting, and 
lossy objects. The geometry of these objects is not easy to set up in a three 
dimensional space. 

In this thesis, the real fractal trees were approximated by planar fractal trees. 
The radar cross section of these trees was calculated by developing a Fortran 
program. The planar fractal trees that were used for this investigation were 
symmetric planar structures composed of perfectly conducting and non—intersecting 


planar dipoles. The geometry of these structures was generated using fractal theory. 


A. CONCLUSIONS 

The radar cross section of the planar fractal trees was calculated using the 
moment method theory. For a given planar structure composed of planar strips and 
illuminated by a plane wave, the program generates the geometry of the PWS 
modes, calculates the impedance between any two of them, and fills the impedance 
matrix {Z]. The voltage on each PWS mode, due to the incident electric field, is also 
calculated, and the voltage column vector [V] is generated. 

The RCS program uses the moment method theory to calculate the current 
distribution on each PWS mode. The radiation equations of electromagnetic theory 
are used to determine the scattered electric field due to the current induced on each 


PWS mode of the structure. The knowledge of both the given incident 
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electric and calculated scattered electric fields leads to the calculation of the radar 
cross section (o/d°), for the given angles 60,y of the incident electric field. 

The calculated radar cross section of a single centered loaded vertical dipole 
was compared with the values given in [Ref. 8, pp. 115], for a similar dipole. 
Although the investigation was for a planar dipole and in [Ref. 8] the dipole was 
cylindrical, the discrepancy of the results was very small. 

In this thesis, five models of planar fractal trees were investigated. The 
geometry of these models was generated by a given program which was modified to 
generate the input data for the developed RCS program. 

The broadside radar cross section, for the monostatic case, was calculated for 
three of these planar fractal trees, for a frequency range of 15 GHz to 75 GHz in 
steps of 15 GHz. This investigation showed that the radar cross section varied in a 
range of 10 dB in the E—Plane. In the H—Plane, the variation of afr was 
symmetric about the 90° axis and smooth at low frequencies and small branching 
angles. For higher branching angles and reduction factors, more variations of o/ Me 
with the frequency were seen in both the E and H planes. 

The maximum RCS at @% = 90°, gy = 0° showed a small variation over a 
frequency range. This range was different on each investigated tree. The variation of 
the maximum RCS followed a straight line at the other frequencies. In each tree and 
at the same frequency, the maximum RCS varied in the range of 3 dB 
approximately. 

It was found that the developed RCS program was not able to calculate the 
radar cross section of fractal trees that were characterized by large values of 
reduction factor and branching angle due to memory restrictions and other 


numerical problems. 
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The scattering from an actual tree will not exactly follow the patterns that 
were described in Chapter 5, but it is felt that the trends would generally remain 


the same. This is specially true for the frequency behavior. 


B. RECOMMENDATIONS 

One recommendation is to investigate the radar cross section of planar fractal 
trees characterized by values other than those used in this thesis. Especially, a 
limited set of reduction factor and branching angle has to be established for the 
same physical length of the initial dipole. 

Another recommendation is to investigate the radar cross section of planar 
fractal trees characterized by the same values of reduction factor and branching 
angles as those that were used in this thesis but with different physical and 
electrical dimensions of the fractal trees. 

In this thesis, the branches of the fractal trees were assumed to be perfectly 
conducting planar strips. The trees were considered without leaves. The radar cross 
section of fractal trees composed of lossy planar strips and with leaves should be 
investigated. 

Finally, the generation of a three dimensional fractal tree and its radar cross 


section may be investigated. 
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APPENDIX A 
PROGRAM RCS 


THIS PROGRAM CALCULATES THE RADAR CROSS SECTION OF A 
PLANAR STRUCTURE COMPOSED OF ARBITRARILY ORIENTED 
PLANAR THIN DIPOLES. THE DIPOLES ARE NOT INTERSECTING. 
THE DIPOLES ARE PERFECT CONDUCTORS. 


GEOMETRY IN Y-Z PLANE 
INPUT DATA: 


F = FREQUENCY IN GHZ 

NW = NUMBER OF DIPOLES 

A.B = COORDINATES OF INCIDENT ELECTRIC FIELD ON Y AND Z 
S25 RESPECTIVELY. 

LW = HALF LENGTH OF EACH DIPOLE IN CM 

WW = HALF WIDTH OF EACH DIPOLE IN CM 

SW, TW = COORDINATES OF THE CENTER OF EACH DIPOLE 

ALONG Y AND Z AXIS RESPECTIVELY 

Porw =‘ANGIE IN. DEGREES BEFWEEN THE Z AXIS (REFERENCE 
AND DIRECTION OF CURRENT FLOW ON EACH DIPOLE 

(POSITIVE ANGLES ARE MEASURED CLOCKWISE) 

NW = NUMBER OF SEGMENTS THAT EACH DIPOLE IS 

SUBDIVIDED 

THITAO, PHIO = ANGLES IN DEGREES OF INCIDENT ELECTRIC 

FIELD. 


OUTPUT DATA: 


NORMALIZED RCS (¢/A2) IN dB 


REAL LS(1:230).S(1:230),T(1:230),PSI(1:230) 

REAL LEN,PG.SU,TU.WID,A,B.MAG(1:230) 

INTEGER I,NMAX,INDX(230),.NT,NP,L,M,N.K,NW 
COMPLEX Z(1:230.1:230).V(1:230),CUR(1:230) 

COMMON/ RC1/ NMAX 

COMMON / RC2/ LS,W1,5,T,PSI 

COMMON/ RC4/ CUR 

OPEN (UNIT=1.FILE='INPUT1',FFORM='FORMATTED') 
OPEN (UNIT=2.FILE='OUTPUT1'. FORM='UNFORMATTED') 
OPEN (UNIT=3.FILE='INPUT' ,FFORM='FORMATTED') 
OPEN (UNIT=4.FILE='OUTPUT'.FORM='FORMATTED') 
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READ(1,*) F,NW,A,B 

PI = 4.* ATAN(1.) 

KO = PI*F/15. 

RAD = PI/180. 

CALL ZMATR (F,NW,Z) 

NP=230 

NT=NMAX 

CALL CLUDCP (Z,NT,NP,INDX) 
READ (3,*,END=11) THITAO,PHIO 
CALL VOLT(F,THITA0,PHIO,A,B,V) 
CALL CLUBSB (Z,NT,NP,INDX,V) 
DO 54 K=1,NT 

CUR(K)=V(K) 

CONTINUE 

CALL RADIAT (F,THITA0,PHIO,NMAX,A,B,RCS) 
WRITE(4,*) THITA0,PHIO,RCS 
PRINT*,RCS 

Core 22 

CLOSE(3) 

STOP 

END 


SUBROUTINE TO COMPUTE THE MUTUAL IMPEDANCE BETWEEN 
THE PWS MODES OF N ARBITRARILY ORIENTED DIPOLES. 


GEOMETRY IN THE Y—Z PLANE 


SUBROUTINE ZMATR Mea. 
REAL F,$(1:230),T(1:230),LG,SG,TG,WI(1:230),PSI(1:230) 
REAL PSIG,PSIW(1:70),LW(1:70),H,W,DH,DW,WIDTH,LENG 
REAL H1,H2,W1,W2,S1,52,T1,T2,PSI1,PSI2,LS(1:230) 
REAL WW(1:70),SW (1:70), TW(1:70),SS(1:20),TS(1:20) 
REAL PG,WID,LEN,SU,TU,A,B 
INTEGER TEMP, P,L,NS(1:70),NMAX,I,N,NW 
INTEGER FUN,G,J,GB,GR,K 
COMPLEX 2Z12,2Z(1:230,1:230) 
COMMON/ RC2/ LS,WI,S,T,PSI 
COMMON/ RC1/ NMAX 
COMMON/ JOHN/ LG,NG,SG,TG,PSIG,LENG 
COMMON/ ZPAR/ H,W,DW,DH 
COMMON/ ZINC/ H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 
DO 1011 =1, NW 

READ(1,*) LW(I), WW(1),SW(1), TW(1),PSIW(I),NS(1) 
CONTINUE 
CLOSE(1) 
NMAX = 0 
DO 75 I=1,NW 

NMAX = NMAX + (NS (I) —1) 
CONTINUE 
P=] 
GB=0 
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DO 50 I=1,NW 

LENG = LW(I) 

LG=2*LW(I)/FLOAT(NS(I)) 

WIDTH=WW(I) 

NG=NS 0 

SG=SWiI 

TG=TW(I) 

PSIG=PSIW(I) 

CALL GEOM(SS,TS) 

GR=NS(I)+GB-1 

M=P 

DO 60 J=1,NS(I)-1 
PSI(M)=PSIG 
LS(M)=LG 
WI(M)=WIDTH 
S(M)=SS(J) 
T(M)=TS(J) 
PRINT*,S(M),T(M) 


M=M+1 
CONTINUE 
GB=GR 
P=GR+41 

CONTINUE 

TEMP=NS(1)-1 

L=1 

G=NMAX 

DO 80 M=1.G 
IF(M.LE.TEMP)GOTO 85 
L=L+1 
TEMP=TEMP+NS(L)-1 
CONTINUE 
K=TEMP 


DO 90 N=M,G 
IF(N.LE.K)THEN 
H=LS(M) 
W=WI(M) 

DW=0 
DH=ABS(M-N)*H 
CALL ZSDIP(F,Z12) 
Z(M,N)=Z12 
PRINT*,Z(M,N) 
WRITE(2) Z(M,N) 
ELSE 

H1=LS(M) 
H2=LS(N) 
W1=WI(M) 
W2=WI(N) 
et 
S2=S(N) 
T1=T(M) 
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T2=T(N) 
PSI1=PSI(M) 
PSI2=PSI(N 
CALL ZPSUR (F,Z12) 
Z(M,N)=Z12 
PRINT*,Z(M,N) 
WRITE(2) Z(M,N) 
ENDIF 
CONTINUE 
CONTINUE 
DO 24 M=2,G 
DO 26 N=1,M-1 
Z(M,N)=Z(N\M) 
WRITE(2) Z(M,N) 
CONTINUE 
CONTINUE 
RETURN 
END 


SUBROUTINE TO COMPUTE THE MUTUAL/SELF IMPEDANCE 
BETWEEN TWO DIPOLES. THE DIPOLES ARE ASSUMED TO BE 
COPLANAR, IDENTICAL AND PARALLEL. THE DIPOLE TO 
DIPOLE IMPEDANCE IS COMPUTED AS THE SUM OF FOUR 
MONOPOLE TO MONOPOLE IMPEDANCES. 


INPUT PARAMETERS: 


F = Frequency of operation (GHz 
H = Half height of the dipole ie 
W = Half width of the dipole (cm 

DH = Longitudinal distance between the two dipoles (cm) 
DW = Transverse distance between the two dipoles (cm) 


SUBROUTINE ZSDIP (F,Z12) 

REAL H, W, DW, DH, F 

COMPLEX 212, ZT 

COMMON/ ZPAR/ H,W,DW,DH 

ZT = (0.,0. 

Z12 = (0.,0.) 

CALL ZSMONP (F, H, W, DW, DH, 0, 1, 0, 1, ZT) 

Z12 = Z12 + ZT 

CALL ZSMONP (F, H, W, DW, DH + H, 0, 1, 1, 0, ZT) 
Tie ioe ae 

CALL ZSMONP (F, H, W, DW, DH —H, 1, 0, 0, 1, ZT) 
fo 2 + ZT 

CALL ZSMONP (F, H, W, DW, DH, 1, 0, 1, 0, ZT) 
Zia 712 27 

RETURN 

END 
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SUBROUTINE TO COMPUTE THE SELF/MUTUAL IMPEDANCE 
BETWEEN TWO IDENTICAL, COPLANAR MONOPOLES. THE 
CURRENT IS ASSUMED TO BE CONSTANT IN THE 
TRANSVERSE DIRECTION. 


REF: R. JANASWAMY, A SIMPLIFIED EXPRESSION FOR THE 
SELF/MUTUAL IMPEDANCE BETWEEN TWO COPLANAR 
AND PARALLEL MONOPOLES, IEEE T—AP, AP-35, 
No. 10, pp. 1174-1176, October 1987. 


INPUT PARAMETERS: 


F = FREQUENCY IN GHz 

H = LENGTH OF EACH MONOPOLE (cm) 

W = WIDTH OF EACH MONOPOLE (cm) 

D = CENTER TO CENTER SPACING BETWEEN THE TWO 
MONOPOLES IN THE DIRECTION TRANSVERSE TO THE 
CURRENT FLOW (cm) 

HH = CENTER TO CENTER SPACING BETWEEN THE TWO 
MONOPOLES IN THE DIRECTION OF CURRENT FLOW (cm) 


111 = TERMINAL CURRENT OF END 1 OF MONOPOLE 1. 
121 = TERMINAL CURRENT OF END 2 OF MONOPOLE I. 
112 = TERMINAL CURRENT OF END 1 OF MONOPOLE 2. 
122 = TERMINAL CURRENT OF END 2 OF MONOPOLE 2. 


NOTE: 111, 121, 112, 122 can assume values only 0 or 1. 

ICODE =0,IF D .LE. 4W 

ICODE = 1, OTHERWISE 

With ICODE = 0, the expression provided in the above paper is used. 
With ICODE = 1, a modified form of the expression provided in the 
above paper is used. (cf. notes) 


OUTPUT PARAMETERS: 


Z12 = COMPLEX IMPEDANCE BETWEEN THE TWO SURFACE 
MONOPOLES. 


SUBROUTINE ZSMONP (F, H 
REAL F, H, W, D, HH, A, B, P 
REAL KW, KH, KD, RC, FR, Fl, 1, 12, 13, 14, SI, Cl 

REAL UABP, ZR, ZI, X, AA (1), BB (1), SQXV, UAB, TINY 
INTEGER 111, 121, 112, 122, M, , NX, KI, ICODE 

COMPLEX 212, AC (—1:1, —1:1), El, E2, Z1, J, CMN, E3, El, FAC 
EXTERNAL FR, FI 

COMMON /PARAM/ N, V, KD, A, B, ICODE 

EI (X) = CI (ABS (X)) — J * SI (X) 

6OxV (x) =SORT (X *X + V* V) 

TINY = 1.E-6 

PI = 4. * ATAN (1.) 

NX =1 


oD, AE i222 212) 
I KO, V, UB, UBP, UA, UAP 


ral 


OAMeieieie® 


KO = PI* F/ 15. 

KW = KO * W 

KH = K0* H 

KD = KO * D 

ICODE = 0 

IF (KD .GT. 4. * KW) ICODE = 1 
KI = 10 


Note that with this choice of KI accurate results are found in the 
region 0 < D < 4 W, and for D > 20 W. In between these two regions, 
accurate results are found with KI = 20. Hence the above choice is 
good only if this code is used for values of D satisfying the above 
inequalities. 


FAC = CMPLX (1., 0.) 
A=KD-2.* KW 

B=KD +2.* KW 

IF (ICODE .EQ. 1) GO TO 2 

AA ( =A 

BB (1) =B 

GOTO: 

AA (1) =-2.* KW 

BB (1) = 2. * KW 

I1 = (121 * COS (KH) — 111) * 112 


2 = (112i COsSthH)) 122 
I3.= "(121 SS COS) 2 
14 = (121 —I11 * COS (KH)) * 122 








J = CMPLX (0.,1.) 

= CEXP (J * KH) 
A@u-1, 1) = 10 4 za 
AC en 


AGC (1, 1) = 13'— 14 aa 
AGG. 1) 18 =a 71 
AC(O,=)) 


= —(AC (211) * 21 4-AC (n=) | Z1) 
AC 0, 1) = (AG Get) * ZAC (si yt 
C= 15. / (2. * SIN (KH) * KW) ** 2 


a9 ao 0) 
DO UM —=—i1 
DOUN——1, 12 


V = Ko * (HH + M* BH) 

IF (ICODE .EQ. 0) GO TO 4 
FAC = CEXP (J*N*V 
CMN = CMPLX (0.,0. a 


SHA 


GOTO5 

UA = SQXV (A) +N*V 
UAP = UA-2.*N*V 
UB = SQXV (B) +N*V 
UBP = UB—2.*N*V 
UAB = SQXV (KD) +N*V 
UABP = UAB—2.*N*V 


El = EI (UB) 
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.GT. TINY) E2 = EI (UA) 


A) 

KD) Gn. TINY) E3 = EI (UAB) 

5*(B*B*El+A*A* E2—2.* KD * KD* E3 + 
CoP (J ~ UB) * (1.4.3 * UBP) + CEXP (—J * UA) * 
Ges J* WAP) —2, *CEXP (— J * UAB) * (1. +11" UABP)) 

CALL HABER (NX, AA, BB, FR, KI, ZR) 

CALL HABER (NX, AA, BB, FI, KI, ZI) 

Z12 = Z12 + A (M, N) * (CEXP (J * N * V) * CMN + (ZR + J * ZI) 

*F 

CONTINUE 

Z12 = Z12 * RC 

RETURN 

END 


REAL FUNCTION FR (XI, NX) 

INTEGER N, NX, CODE 

REAL XeV, Tl, KD, ApS, XINX), TKW)T2, CI, SI, TINY 
COMPLEX EI, J 

COMMON ans ix) N, V, KD, A, B, CODE 
EI (X) = CI (ABS (X)) — J * SI (X) 

TINY = 1.E-6 

X = XI (1) 

J = CMPLX (0.,1.) 

iW (CODE .EQ.1)GOTO1 

Mm — SORT (X* X+ V* V) 

FR = COS (T1) 

ABS (V).GT. TINY) FR =FR*(1.—N* V/T1) 
teen CE hy bR = FR * A 

(Xx GT. KD) PR =—FR* B 

GO TO 2 

TKW = 0.5 * (B— A) 

2 — SORT AD + X) * (KD + X) + V * V) 
FR = REAL (El (T2 + N* YY) 

FR = FR * (TKW — ABS (X) 

END 





REAL FUNCTION FI (XI, NX) 

INTEGER N, NX, CODE 

REAL X, V, Tl, KD, A, B, XI (NX), TKW, T2, CI, SI, TINY 
COMPLEX EI, J 

COMMON /PARAM N, V, KD, A, B, CODE 
EI (X) = CI (ABS (X)) —J * SI (X) 

X = XI (1) 

TINY = 1.E 

J = CMPLX (0., 1.) 

Tl = SQRT (X*X+V*V) 

IF (CODE .EQ. 1)GOTO1 
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FI = —SIN (T1) 

IF (ABS (V) .GT. Ne FI = FI* (1.-N* V/T1) 
IF (X .LE. KD) FI=FI* A 

10 (OC GHES 10) el ae =B 

GO 102 

TKW = 0.5* (B—A) 

T2 = SORT ((KD + X) * (KD + X)4 VV) 

FI = AIMAG (EI (T2 + N * V)) 

FI = FI * to ABS (X)) 

END 


SUBROUTINE TO COMPUTE THE COORDINATES OF THE PWS 
MODES OF EACH DIPOLE IN THE Y-Z PLANE 


SUBROUTINE GEOM (S,T) 

REAL THI,H1,H2,Y,Z,SG,TG,LG,PSIG 

REAL P,Q,S(1:20),T (1:20), LENG 

INTEGER M,NG.K 

COMMON/ JOHN/ LG,NG,SG,TG,PSIG,LENG 

COSD(X) = COS(X*RAD) 

SIND(X) = SIN(X*RAD) 

PIl=4.*ATAN(L. 

RAD=PI/180. 

THI=90.—PSIG 

H1=COSD(THI) 

H2=SIND(THI) 

Y=SG-(LENG*H1 

Z=TG-(LENG* H2 

P=LG*H1 

Q=LG*H2 

S(1)=Y+P 

T(1)=Z+Q 

DO 225 K=2,NG-1 
S(K)=S(1)+(K—1)*P 
T(K)=T(1)+(K-1)*Q 

CONTINUE 

RETURN 

END 





SUBROUTINE ZPSUR (F,Z12) 

REAL H1, W1, H2, W2, PSI, F, YSTAR, ZSTAR, KO, AI (2), BI (2) 
REAL C (3), D (3), RT, ZT, COT, CSEC, KOS, X, PI, PSI1, PSI2 
REAL S1, T1, $2, T2, RINTG, IMINTG, RESULT, SIND, COSD,SI,CI 
INTEGER NX, KI 

COMPLEX Z12, J 

EXTERNAL RINTG, IMINTG 

COMMON/ ZINC/ H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 

COMMON /PARAM2/ C, D, RT, ZT, CSEC, COT, KOS, J, KO 
SIND (X) = on Co - / 180.) 
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COSD (X) = COS (X * PI / 180.) 

DATA Al, BI /2* Se 

PI = 4. * ATAN (1.) 

PSI = PSI2 — PSII 

IF (ABS (PSI).LE.0.4) PSI= SIGN(1.,PSI)*0.4 
J = CMPLX (0.,1.) 

NX =2 
KI = 3 
K0 = PL : 


180 
/ 


F 
{gn (KO * H1) 


SS ® cos (Ko" H1) * C (1) 
1. / SIN (KO * H2) 
D (1) 
—2. * COS (K0'#/H2) *'D (1) 
CSEC = 1. / SIND (PSI) 
KOS = COSD (PSI) 
COT = KOS * CSEC 
YSTAR = ($2-S1) * COSD (PSI1) — (T2-T1) * SIND Ea 
ZSTAR = (S2-S1) * SIND (PSI1) + (T2-T1) * COSD (PSI1 
RT = YSTAR * CSEC — H? 
ZT = YSTAR * COT —H1-ZSTAR 
CALL HABER (NX, Al, BI, RINTG, KI, RESULT) 
Z12 = RESULT 
CALL HABER (NX, Al, BI, IMINTG, KI, RESULT) 
Wi 7124 J? RESULT 
Z12 = -3.75 * Z12 
RETURN 
END 


Oe 2oO0@ 
Nw wwe 


REAL FUNCTION SI (XI) 

REAL XI 

REAL AF (4), BF (4), AG (4), BG (4), X, X2, X4, X6, X8, FX, GX, 
1, SGN 


DATA AE / 38.027264, 265.187033, 335.67732, 38.102495 / 
DATA BF / 40.021433, 322.624911, 570.23628, 157.105423 / 
DATA AG / 42.242855, 302.757865, 352.018498, 21.821899 / 
DATA BG / 48.196927, 482.485984, 1114.978885, 449.690326 / 
SI = 0. 

IF (XI .EQ. 0.) RETURN 

SGN = +1. 

[P(e LT. 0.) S@n= — 
X = ABS (XI) 
x2 =x 


X6 = X4 * X2 
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X8 = X4* X4 

IF (X .LT. 1.) THEN 

SI = X*(1.— X2/ 18. + X4 / 600. —X6 / 35280. + X8 / 3265920.) 

SI = SI * SGN 

ELSE 

FX = (1. + AF (1) / X2 + AP @) / x4 AF @)/ Xo 
X8 


FX = PO O* (1 + BF (1) / X2 + BF (2) / X4 + BF (3) / X6+ 


BF 
ce (FAG (1) X2 + AG (2) / X44 AG (8) / X64 AG (4) 
Gx Ox (x2 2* (1. + BG (1) / X2 + BG (2) / X4 + BG (3) / X6+ 
BG(4) / X8)) 
PI = 4. * ATAN (1. 


= x= AINE] 0+ P * Pl)) * 2. * PI 

SI = SGN * (PI / 2.—FX * Gos (x) — GX * SIN (X)) 
END IF 

END 


REAL EP BE { a ae 

REAL AF (4), B 

REAL PI 

DATA AF / 38.027264, 265.187033, 335.67732, 38.102495 / 

DATA BF / 40.021433, 322.624911, 570.23628, 157.105423 / 

DATA AG / 42. 242855, 302. 757865, 352. 018498, 21.821899 / 

DATA BG / 48.196927, 482.485984, 1114. 978885, 449.690326 / 

IF (X .LE. 0.) THEN 

PRINT *, ‘Invalid argument for CI (x)','x =', x 

RETURN 

ELSE 

x2 aes 

NS 2 Xe 

IF (X .GE. 20.) THEN 

FX = 1 x (ee, 

GX =1. / X2* (1.-6. / X2) 

Goon 

END IF 

XG = X4* X2 

X8 = X4 * X4 

IF (X .LT. 1.) THEN 

CI = 0.57721566 + ALOG (X) — X2 * (0.25 — X2 / 96. + X4 / 4320. 
— X6 / 322560.) 

ELSE 


PX = (1. + AF (1) / X2 + AF (2) / X4 + AF (8) / X6 + AF (4) / 

8 

FX = FX / (X * (1. + BF (1) / X2 + BF (2) / X4+ BF (3) / X6+ 
Be 

GX = (1. + AG ( 1) /X2+ AG (2) / X44 AG @)/ X6 + AG) 


GX = GX / (X2 * (1. + BG (1) / X2 + BG (2) / X4 + BG (3) / X6+ 


); BG Gy G2, KG eee 
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BG(4) / X8)) 

PI = 4.7 AT oe) 

M2 Pit 2" el 
G XCOS (X 


X = X-AINT 

Cl = FX * SIN (X)— ) 
END IF 

END IF 

END 


REAL FUNCTION RINTG (X, NX) 
COMPLEX J, El, TERM1, TERM2, INTGR 
INTEGER NX, M, N, P, Q 
REAL * 4 PSil, PSI2, $1, T1, $2, T2 
RUT * 44C (3), D (3), R, Z, RT, 20, U, V, CSEC, COT, KOS, Y 
REAL * 4 X (NX), Hl, H2, W1, W2, KO, PZQR, SI, CI 

,TT, ZM, RN, RMN, TESC, TESD, PI, PZQRP 
COMMON/ ZINC/ H1,H2,W1,W2,S1,82,T1,T2,PSI1,PSI2 
COMMON /PARAM2/ C, D, RT, ZT, CSEC, COT, KOS, J, KO 
EI (Y) = Cl (ABS(Y)) — J * SI (Y) 
Y= 


DO3M=1,3 
IF (MEQ, 2 AND. TESC .LE. 1.E-6) GO TO 3 
= (M-1) * H1 
IM SU» Wi OO Vi Wo CSEC + ZT +7, 
TERM2 = (0.,0.) 
DO2N=1,3 
IF (N .EQ. 2 .AND. TESD .LE. 1.E-6) GO TO 2 
R = (N-1) * H2 
RN =-U * Wi * CSEC + V* W2* COT+RT+R 
IF (ABS (KO * RN) .LE. 0.5E-1) THEN 
PZQRP = KO * ABS (ZM) — AINT (KO * ABS (ZM) / (2. * PI)) * 


2." PI 
TERM1 = CEXP (—J * PZQRP) * TT 
GO TO 4 
END IF 


IF (ABS (KO * ZM) .LE. .0.5E=1) THEN 
PZQRP = KO * ABS (RN) — AINT (KO * ABS (RN) / (2. * PI) * 
* Pj 


TERM) = CEXP (=J * PZORP) * TT 

GO TO4 

END IF 

RMN = SQRT (RN * RN + ZM * ZM —2. * RN * ZM * KOS) 
TERM1 = (0.,0.) 

D@1 P =-1, 1,2 

DOT O=—1, 1,2 

PZQR =P*ZM+Q*RN 
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PZQRP = PZQR * KO — AINT (PZQR * KO / (2. * PI)) * 2. * PI 
TERMI = TERM1 + P * Q * CEXP (J * PZQRP) * 
EI (KO * (RMN + PZQR)) 
CONTINUE 
TERM2 = TERM2 + D (N) * TERMI 
CONTINUE 
INTGR = INTGR + C (M) * TERM2 
CONTINUE 
RINTG = REAL (INTGR) 
END 


REAL FUNCTION IMINTG (X, NX) 
COMPLEX J, El, TERM1, TERM2, INTGR 
INTEGER NX, M, N, P,Q 
REAL * 4 PSI1, PSI2, $1, T1, $2, T2 
REAL * 4 C (3). DG), R, 2, RI ZU, V. CSEO, COM KOsmn 
REAL * 4 X (NX), H1, H2, W1, W2, KO, PZQR, SI, CI 

TT, ZM, RN, RMN, TESC, TESD, PI, PZQRP 
COMMON/ ZINC/ H1,H2,W1,W2,81,$2,T1,T2,PSI1,PSI2 
COMMON /PARAM2/ C, D, RT, ZT, CSEC, COT, KOS, J, KO 
EI (Y) = Cl (ABS(Y)) —J * SI (Y) 
U=X(l 
V=X(2 


INTGR = (0.,0.) 


)/(.—KOS))) 


DO3M=1,3 

IF (M .EQ. 2 .AND. TESC .LE. 1.E-6) GO TO 3 

Z = (M-1) * Hl 

7M =—U * W1 * COT + ¥* W2* CSEC par + Z 

TERM2 = (0.,0.) 

DO2N=1,3 

IF (N .EQ. 2.AND. TESD .LE. 1.E-6) GO TO 2 

R = (N-1) * H2 

RN = —U * WI * CSEC + Vi" Woe Colne 

IF (ABS (KO * RN) .LE. 0.5E-1) THEN 

PZQRP = KO * ABS (ZM) — AINT (KO * ABS (ZM) / (2. * PI)) * 
ey 


TERM1 = CEXP (-J * PZQRP) * TT 

GO TO 4 

END IF 

IF (ABS (KO * ZM) .LE. .5E-1) THEN 

PZQRP = KO * ABS (RN) — AINT (KO * ABS (RN) / (2. * PI)) * 
2.* PI 

TERM1 = CEXP (-J * PZQRP) * TT 

GO TO 4 

END IF 

RMN = SQRT (RN * RN + ZM * ZM —2. * RN * ZM * KOS) 
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TERMI = (0.,0.) 

DO MP’ = -1, 1, 2 

DO1Q=-1,1,2 

PZQR = P*27M+Q*RN 

PZQRP = Ko * PZQR — AINT (KO * PZQR / (2. * PI)) * 2. * PI 

TERMI = TERMI + P * Q * CEXP (J * PZQRP) * 
EI (KO * (RMN + PZQR)) 

CONTINUE 

TERM2 = TERM2 + D (N) * TERMI 

CONTINUE 

INTGR = INTGR + C (M) * TERM2 

CONTINUE 

IMINTG = AIMAG (INTGR) 

END 


Subroutine to compute a sequence of estimates EST1 (K) and 

EST2 (K), 1 .LE. K1 .LE. K .LE. K2 for the N—dimensional integral 
Bl BN 
Weteemlitte FUN (x1, x2... XN) dx! dx2.., dx3 
Al AN 

by Haber's method. 

Ref: P.J. Davis and P. Rabinowitz, Methods of Numerical 

Integratation, Academic Press, 1984. 


For each estimate EST1 (K), two additional quantities ERR1(K) 

and DEV1 (K) are computed. If the values of DEV1 pi do not 

vary by more than 10% between consecutive values of K, then 

ERR1 (K) can be taken as a reliable bound on the difference 

between EST1 and the integral. A similar situation holds for 

EST2, DEV2, and ERR2 (K). The total number of functional 
evaluations is 4 * (K1 ** N + (K1+1) ** N+... + K2 ** N) and K2 
should be chosen so as to make this may be halfed by eliminating the 
computation of the EST2 (K). In other situations, these values 

are much better than the EST1 (K). A program FUNCTION FUN (X, N) 
must be supplied by the user with X declared by the statement 
DIMENSION X (N). FUN must be declared EXTERNAL in the calling 
program. If N<1lorN > 10or Kl <1 or K2 < K1, the program 
terminates with IND = 0. Otherwise IND = 1 


Modified by R. Janaswamy so that the output is average of EST1 
EST2. Also, K1 = K2 = K. 


SUBROUTINE HABER (N, LL, UL, FUN, K, RESULT) 
INTEGER N, IND, KEY, I, K, J 
DOUBLE PRECISION AL (10), BE (10), GA (10), B, G 
REAL FUN, Y1, Y2, Y3, Y4, EST1, EST2, ERR1, ERR2, 
DEV1, DEV2, RESULT 
REAL * 8 Si, $2, Dl, D2 
REAL LL (N), UL (N), DEX (10), P1 (10), P2 (10), P3 (10), P4 (10), 


& QI (10), Q2 (10), Q3 (10), Q4 (10), RAN (10), AKN, AK, T, JAC 


REAL AKI, BK 


iS 


Sa eS 


EXTERNAL FUN 

DATA AL/.4142135623730950, .7320508075688773, .2360679774997897, 
6457513110645906, .3166247903553998, .6055512754639893, 
.1231056256176605, .3589989435406736, .7958315233127195, 
3851648071345040/ 

IND = 0 

IF (N .LT. 1 .OR. N .GT. 10) RETURN 

IND = 1 


AKN = AK ** N 

T = SQRT (AKN) * AK 

BK = 1. / AK 

KEY = KEY +1 

IF (KEY .EQ. 1) GO TO 6 

KEY = KEY - 1 

J=1 

IF (DEX (J)LeT kineo 105 
DEX (J) = DEX (J) + 1. 

GO TO 6 


=] 


Qe eS eee 


P4 (I Is GbE) 1.— GA (1) - BK 
Q4 (I) = LL (I) + RAN (1) * P4 (1 

Y1 = FUN (QI, N 

Y2 = FUN (Q2, N 

Y3 = FUN (Q3, N 

Y4 = FUN (Q4, N 


Sec) + Yl + Y2 

bi = Di + (Y1 — Y2) ** 2 

S2 =S2+ Y3+ Y4 

ie — 2 + (Vl + Y3 — Y2— Y4) 2 
GO TO5 

EST] = 0.5 * $1 / AKN 

ERRI = 1.5 * DSQRT (D1) / AKN 
DEV1 = ERR] * 

EST2 = 0.25 * G1 + $2) / AKN 
ERR2 = 0.75 * DSQRT (D2) / AKN 
DEV2 = ERR2* T* AK 

RESULT = 0.5 * (EST1 + EST2) * ABS (JAC) 
RETURN 

END 


SUBROUTINE VOLT WHICH CALCULATES THE VOLTAGE MATRIX 
VM 


GEOMETRY IN THE Y—Z PLANE 


SUBROUTINE VOLT(F,THITAO,PHIO,A,B,V) 
REAL TSI,TCO,PS,PC,KY,KZ,A,B,KO0,F,KX 
REAL MAR,GTI,PRS,PRC,MOD,KZHTA 
REAL S(1:230),T(1:230) ,PS1(1:230),LS(1:230),W1(1:230), THITA0,PHIO 
INTEGER NMAX 
COMPLEX V(1:230),J,STR,VOM 
COMMON/ BRAVO/ KX,KY,KZ 
COMMON/ RC2/ LS,WLS,T,PSI 
COMMON/ RC1/ NMAX 
COSD(X)=COS(X*RAD) 
SIND(X)=SIN(X*RAD) 
PI=4.*ATAN(I.) 
RAD=PI1/180. 
K0O=PI*F/15. 
J=CMPLX(0.,1.) 
TSI=SIND(THITAO 
TCO=COSD(THITA0) 
PS=SIND(PHIO) 
PC=COSD(PHIO0) 
KY=KO*TSI*PS 
KX=KO0*TSI*PC 
KZ=K0*TCO 
DO 234 K=1,NMAX 

PRS=SIND(PSI(K)) 
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PRC=COSD(PSI(K)) 

MAR=KY*S(K)+KZ*T(K) 

STR=CEXP(—J*MAR) 

VOM=(STR/SIN(LS(K)*K0))*(A*PRS+B*PRC) 

KZHTA=KY*PRS+KZ*PRC 

MOD=COS(K0*LS(K))—COS(K ZHTA*LS(K)) 

GTI=KZHTA**2-K0**2 

IF(ABS(K ZHTA-K0)/KO.LT.1.E-2.0R. 
ABS(KZHTA+K0)/K0.LT.1.E—-2) THEN 

V(K)=VOM*LS(K)*SIN(KO*LS(K)) 

ELSE 

V(K)=(2*K0/GTI)*VOM*MOD 

ENDIF 


CONTINUE 
RETURN 
END 


SUBROUTINE CLUDCP (A, N, NP, INDX) 
INTEGER NP, N, INDX (NP), I, J, K, P, NMAX 
PARAMETER (NMAX = 230) 
COMPLEX A (NP, NP), TEMP, ETA, W (NMAX) 
REAL AAMAX, DUM 

DO 1K =1, N-1 

AAMAX = CABS (A (K, K)) 
P=K 

DO 2I=K+41,N 

DUM = CABS (A (I, K)) 

IF (DUM .GT. AAMAX) THEN 
AAMAX = DUM 

P=I 

END IF 

CONTINUE 

INDX (K) = P 

DO3J=1,N 

TEMP = A (K, J) 

A (K, J) =A (P, J) 

A (P, J) = TEMP 

CONTINUE 

DO 4J=K+1,N 

W (J) =A (K, J) 

CONTINUE 

DO 51=K+1,N 

ETA = A (I, K) / A (K, K) 

A (I, K) = ETA 
DO6J=K+1,N 

A (1, J)=A CL, 3) EPA) 
CONTINUE 

CONTINUE 

CONTINUE 

RETURN 

END 
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SUBROUTINE CLUBSB (A, N, NP, INDX, X) 
INTEGER N, NP, I, J, INDX (NP), NMAX 
PARAMETER (NMAX = 230) 

COMPLEX A (NP, NP), X (NP), TEMP, Y (NMAX) 


DO PERMUTATIONS ON THE EXCITATION VECTOR USING THE 
INFORMATION ON THE ROW OPERATIONS DONE IN CLUDCP. 


DO1K=1,N-1 
TEMP = X (K) 

X (K) = X (INDX (K)) 
X (INDX (K)) = TEMP 
CONTINUE 


FORWARD ELIMINATION 


BACK SUBSTITUTION 


DO31=N,1,-1 

X (I =Y (1) 
DO4J=I+1,N 

X ().= 240) — A (1, JX (J) 
CONTINUE 

X (I) =X (1) / A (1,1) 
CONTINUE 

RETURN 

END 


SUBROUTINE RADIAT TO COMPUTE THE RADAR CROSS SECTION 
OF A PLANAR STRUCTURE COMPOSED OF ARBITRARILY 
ORIENTED PLANAR DIPOLES. 


SUBROUTINE RADIAT (F,THITA0,PHI0,NMAX,A,B,RCS) 
REAL RC,THITAO,PHIO,LS(1:230), W1(1:230) ,S(1:230),T(1:230),P 
REAL CR,F,KO0,PI,UN,DM,EM,RAD,P,R,A,B, THITA,PHI,RCS PSI(I: 230) 
REAL KX,KY,KZ,KRON,KG 

INTEGER G,NMAX 

COMPLEX J,NTHIO,NPHIO,TEMP1,TEMP2,CUR(1:230) 
COMMON/ BRAVO/ KX,KY,KZ 

COMMON/ RC2/ LS,W1,8,T,PSI 

COMMON/ RC4/ CUR 

COSD(X) = COS(X*RAD) 

SIND(X) = SIN(X*RAD) 
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PI = 4. * ATAN(1.) 

RAD = PI/180. 

KO = PI*F/15. 

UN=120.*PI 

THITA =180.— THITAO 

PHI = 180. + PHIO 

J=CMPLX(0.,1.) 

TE 

TEMP2=(0.,0. 

NTHIO=(0.,0. 

NPHIO=(0.,0. 

DO 58 G=1,NMAX 
P1=PSI(G) 
DM=K0*(S(G)*SIND(THITA)*SIND(PHI)+T(G)*COSD(THITA)) 
EM=K0*(SIND(P1)*SIND(THITA)*SIND(PHI)+ 

+COSD(P1)*SIND(THITA)) 
P=SIND(P1)*COSD(THITA)*SIND(PHI)—COSD(P1)*SIND(THITA) 
TEMP1=P*CEXP(J*DM)*CUR(G)/SIN(K0*LS(G)) 
TEMP2=CEXP(J*DM)*CUR(G)*SIND(P1)*COSD(PHI)/SIN(K0*LS(G)) 
nae /K0.LT.1.E—2.0R.ABS(EM+K0)/K0.LT.1.E—2) THEN 
KGELS G)*SIN(K0*LS(G)) 
ELSE 


KG=—2*K0*(COS(EM*LS(G))—COS(K0*LS(G)))/(EM**2—K0**2) 
ENDIF 
NTHIO=TEMP1*KG+NTHIO 
NPHIO=NPHI0+KG*TEMP2 
CONTINUE 
CR=CABS(NTHIO)**2+CABS(NPHIO)**2 
KRON=(A KY 4B°KZ)/KX 
RC=(K0**2)*(UN**2)*CR/(4.*PI*(ABS(A)**2+ABS(B)**2+ 
+ABS(KRON)**2)) 
IF((RC/((30./F)**2)).LT.1.E-6) THEN 
RCS = 10.*ALOG10(1.E-6) 
ELSE 


RCS=10.*ALOG10(RC/((30./F)**2)) 
ENDIF 
RETURN 
END 
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APPENDIX B 
PROGRAM TREE 


PROGRAM TREE 


PROGRAM TO GENERATE A PLANAR FRACTAL TREE WITH 
REDUCTION FACTOR (COND) AND BRANCHING ANGLE THETA 
WRITTEN BY T.R. NELSON, PhD, UNIVERSITY OF CALIFORNIA 
SAN DIEGO, LA JOLLA, CA, 90293, AND MODIFIED BY 

LT. JOHN DEMIRIS TO GENERATE THE INPUT DATA FOR 

THE RCS PROGRAM. 


INPUT DATA: 


L = NUMBER OF BRANCHING LEVELS. 
N = NUMBER OF BRANCH SEGMENTS. 
ILEN = INITIAL LENGTH. 

ISP = INITIAL STARTING POINT. 
COND = REDUCTION FACTOR. 
THETA = BRANCHING ANGLE. 


OUTPUT DATA: 


IX; TY = COORDINATES OF FIRST AND END POINT OF EACH 
GENERATED BRANCH. 
INPUT DATA FOR RCS PROGRAM 


REAL ZZ (1:5,1:1024) ISP ,ILEN,IX,lY ITH1 

OPEN (UNIT=1, FILE = ‘INGEOM',FORM = 'FORMATTED') 
OPEN (UNIT=2, FILE ="OUT', FORM = ‘'FORMATTED') 
CEE ONIT=10, FILE =INPUTI',FORM='FORMATTED'") 
READ(1,*,END=10) L,N,ILEN,ISP,COND,THETA 





Z2N=N 

PI = 3.14159265 
IX = 0. 

IY = 0.-ILEN 
IX = zLAe, | 
IY = Z2(1,1 

DO 100 LOOP =1,L 
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ZLOOP = LOOP 

NPTS = N**ZLOOP 

INC =0 

NPREV = N**(ZLOOP-1) 

PREV = NPREV 

ANGLE = (ZN-1.0)/2.0 

ZINDX = -ANGLE 

DO 200 J =1,NPTS 

ZJ=J 

IF (INC.GE.N) ZINDX=—-ANGLE 
IF (INC.GE.N) INC=0 
IW=NPTS—J+1 

IR = PREV-ZJ/ZN+1.0 
ZL=ZZ(4,1R) 

T = Z2Z(3,IR 
X = ZZ(1,IR 
Y = Z22(2,1R 
ZL1 = ZL*COND 

IF (ZL1.LT.0.1) GOTO 300 
TH1 = ZL1/33.0 

Tl = T+THETA*ZINDX 
ZZ(3,1W) = Tl 

T2 = T1*PI/180.0 

IX=Y 

IY =X 

TEMPX1=IX 
TEMPY1=lY 


IX, TY ARE THE COORDINATES OF THE FIRST POINT OF THE 
BRANCH ALONG THE X,Y AXIS RESPECTIVELY 








WRITE(2,*) IX,IY 
X= ZL cos) x 
Y1 = ZLI*SIN(T2)+Y 
IX=Y1 

1Y=X1 


IX,TY ARE THE COORDINATES OF THE END POINT OF THE 
BRANCH ALONG THE X,Y AXIS RESPECTIVELY 


WRITE(2,*) IX,IY 
TEMPX2=1X 
TEMPY2=lY 


S AND T ARE THE COORDINATES OF THE CENTER OF EACH 
BRANCH 


S = (TEMPX2+TEMPX1)/2. 

T = (TEMPY2+TEMPY1)/2. 
PSI1=ATAN((TEMPX2—TEMPX1)/(TEMPY2-TEMPY1)) 
PSI=PSI1*180. /PI 
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ZLHALF=ZL1/2. 
THIHALF=TH1/2 
WRITE(10,*) ZLHALF,THIHALF,S,T,PSI 
ZZ(1,1W)=X1 
ZZ(2,1W)=Y1 
ZZ(4.1W)=ZL1 
ZINDX=ZINDX+41.0 
INC=INC+41 

200 CONTINUE 

100 CONTINUE 

300 STOP 

END 
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